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ABSTRACT 


A  numerical  algorithm  suitable  for  rapid  numerical  solution  of 
advanced  turbulence-model  equations  on  elliptic  regions  has  been 
devised.  The  algorithm  is  particularly  useful  for  obtaining  accu¬ 
rate  solutions  when  (a)  one  of  the  dependent  variables  is  singular 
approaching  a  solid  boundary  and  (b)  the  number  of  mesh  points 
available  to  resolve  the  flowfield  is  small.  The  new  algorithm 
is  used  in  conjunction  with  the  MacCormack  semi-implicit  scheme 
for  solving  the  time-averaged  Navier-Stokes  equations  to  obtain 
solutions  for  a  laminar  Mach  2  boundary  layer  and  for  a  turbulent 
Mach  3  boundary  layer.  In  both  cases,  numerical  accuracy  is  shown 
to  be  excellent  for  meshes  having  less  than  20  points  normal  to 
the  surface.  As  a  corollary  result  of  this  study  we  have  shown 
that,  using  the  new  algorithm  and  the  MacCormack  scheme,  time- 
averaged  Navier-Stokes  solutions  are  now  feasible  using  a  micro¬ 
computer. 


FOREWORD 


This  report  presents  a  summary  of  research  performed  in  Contract 
F49620-78-C-0024  during  the  period  March  1,  1978  through  December 
31,  1981.  This  research  was  sponsored  by  the  Air  Force  Office  of 
Research  (AFSC),  United  States  Air  Force.  The  Air  Force  program 
monitor  was  Dr.  James  Wilson. 

Study  participants  were  Dr.  David  C.  Wilcox,  principal  investi¬ 
gator,  and  Barbara  A.  Wilcox,  contract  administrator  and  data 
processing  support.  Manuscript  preparation  was  done  by  Kinley  A. 
Wilcox. 
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1.  INTRODUCTION 

During  the  course  of  this  research  project  (and  also  in  a  numerical 
study  funded  by  the  NASA  Ames  Research  Center"'")  we  have  found  that 
conventional  time-marching  numerical  methods  encounter  serious 
difficulty  in  solving  advanced  turbulence-model  equations.  This 
j  difficulty  has  impeded  our  progress  toward  the  ultimate  goal  of  this 

research  project,  viz,  development  of  a  numerical  method  for  simu¬ 
lating  viscous  flow  in  axial  turbomachines.  Thus,  during  the  past 
year  our  efforts  have  focused  on  the  cause  of  and  the  resolution  of 
|  this  numerical  difficulty. 

Prior  to  this  study  we  found  that  numerically-accurate  steady-flow 
solutions  generally  cannot  be  obtained  for  complex  turbulent  flows 
under  the  following  conditions: 

^  1.  an  advanced  turbulence  model  with  additional 

dependent  variables,  one  of  which  usually  is 
singular  near  solid  boundaries,  is  used; 

2.  a  severely  limited  number  of  points  are  avail- 

i  able  to  define  the  computational  mesh  for  the 

1  region  of  interest; 

3.  the  equations  are  integrated  through  the  viscous 
sublayer. 

i  For  our  axial  turbomachine  applications,  all  three  of  the  conditions 
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cited  above  exist.  Firstly,  we  are  using  the  Wilcox-Rubesin  two- 

equation  model  of  turbulence  which  involves  a  dependent  variable 
that  becomes  singular  upon  approach  to  a  solid  boundary.  Secondly, 

»  the  flowfield  in  axial  turbomachines  inherently  is  three  dimensional 

ft 

and,  because  computers  have  limited  memory,  we  anticipate  having  no 
more  than  15  mesh  points  normal  to  a  blade  surface  contained  within 
boundary  layers.  Finally,  because  we  feel  that  physically  meaning- 
v  ful  solutions  on  separation  regions  require  integration  through  the 

viscous  sublayer,  and  because  we  want  our  method  to  treat  the  sepa¬ 
rated  case,  our  applications  include  integration  through  the  viscous 
sublayers. 

Prior  to  this  study,  Wilcox^"  identified  two  possible  sources  of 
inaccuracy  and/or  numerical  Instability  which  might  plague  conven¬ 
tional  time-marching  methods.  He  demonstrated  that  the  singular 
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behavior  of  the  turbulent  dissipation  rate  near  a  solid  boundary 
cannot  be  accurately  computed  with  central  differences.  In  an 
unsuccessful  attempt  at  analytically  removing  the  singularity  with 
a  straightforward  change  of  dependent  variables,  improved  numerical 
accuracy  was  accompanied  by  slowly  growing  numerical  oscillations 
elsewhere  in  the  flowfield.  Wilcox  speculated  that  the  source  of 
these  oscillations  was  the  numerical-algorithm-induced  (and  physi¬ 
cally  unrealistic)  coupling  between  the  turbulence  production 
processes  and  the  unsteady  rate  of  change  of  the  turbulence. 

In  this  study,  we  have  devised  a  numerical  algorithm  which  lacks  the 
shortcomings  of  conventional,  central-differences-oriented,  time¬ 
marching  methods  when  applied  to  turbulence-model  equations.  The 
method's  accuracy  and  stability  have  been  assessed  in  two  viscous- 
flow  computations  (both  performed  on  a  microcomputer),  neither  of 
which  can  be  done  accurately  with  a  conventional  time-marching 
method  unless  far  more  mesh  points  are  used. 

Section  2  presents  the  equations  of  motion  and  pinpoints  the  root 
causes  of  the  numerical  inaccuracy  and  instabilities  encountered  in 
our  attempts  to  compute  viscous  flowfields.  In  Section  3S  a  new 
algorithm  is  devised  to  circumvent  these  difficulties;  application  to 
a  simplified  turbulence-model  equation  demonstrates  a  dramatic  impro¬ 
vement  in  numerical  accuracy  over  conventional  algorithms.  Section 
4  presents  results  of  two  time-averaged  Navier-Stokes  computations, 
one  for  a  laminar  Mach  2  boundary  layer  and  the  other  for  a  turbulent 
Mach  3  boundary  layer.  In  both  cases,  results  are  (correctly) 
obtained  which  appear  unobtainable  with  a  conventional  time-marching 
method.  The  concluding  section  summarizes  results  and  conclusions. 


The  Bibliography  and  Contract  Overview  sections  summerize  publications 
and  highlights  of  the  overall  research  effort  in  this  Contract. 


2.  EQUATIONS  OF  MOTION  AND  ASSOCIATED  NUMERICAL  DIFFICULTIES 


This  section  first  presents  the  basic  equations  of  motion  used  in 
this  project.  A  summary  of  the  numerical  difficulties  encountered 
in  trying  to  solve  these  equations  follows,  including  a  discussion 
of  truncation  error  and  numerical  stability. 

2.1  EQUATIONS  OF  MOTION 

The  equations  of  motion  used  in  all  of  our  computations  are  those 

p 

devised  by  Wilcox  and  Rubesin  .  The  model  assumes  the  Reynolds 
stress  tensor,  Tjjj  is  proportional  to  the  mean  strain-rate  tensor, 
S^j ,  as  follows. 

Tij  =  2pe  (Sij  "  I  3x^  6ij}  "  f  pe  6ij  (1) 


In  equation  (1),  p  is  mass  density,  e  is  the  eddy  diffusivity,  u^  is 
the  velocity  vector,  is  position  vector,  5^  is  the  Kronecker 
delta,  and  e  is  the  turbulent  mixing  energy.  All  quantities  are 
understood  to  be  mass  averaged  in  the  sense  introduced  by  Favre  . 

The  mean  conservation  equations  are: 
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where  t  denotes  time,  p  is  static  pressure,  y  is  molecular  viscosity, 
k  is  the  sum  of  molecular  and  turbulent  heat  dif fusivities ,  T  is 
static  temperature,  o#  is  a  turbulent  closure  coefficient  which  will 


be  discussed  below,  and  E  is  the  total  energy  defined  by  the  follow¬ 
ing  relation. 

E  =  CyT  +  +  pe  (5) 


where  Cy  is  specific  heat  coefficient  at  constant  volume. 


In  addition  to  Equations  (1-5),  the  Wilcox-Rubesin  model  introduces 
two  additional  "rate  equations"  for  the  turbulent  mixing  energy,  e, 
and  the  turbulent  dissipation  rate,  <u,  from  which  the  eddy  diffusivity 
is  computed.  The  equations  are  as  follows. 
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The  eddy  diffusivity,  e,  and  the  turbulent  length  scale,  are 
obtained  from  e  and  oj  according  to: 


e  =  y*  e/w 
l  =  e^/u 


(8) 

(9) 


Finally,  we  introduce  a  turbulent  Prandtl  number,  PrT,  and  express  the 
turbulent  heat  flux  vector  as  being  proportional  to  the  mean  temperature 
gradient  so  that  the  thermal  conductivity,  k,  in  Equation  (5)  becomes: 

k  =  }  P/PrL  +  pe/PrT  }  Cy  (10) 


Where  PrT  is  laminar  Prandtl  number.  Six  closure  coefficients  appear 

U 

in  Equations  (6-8),  viz,  6,  3* ,  y,  Y*,  a,  and  a*.  The  values  of  these 
closure  coefficients  have  been  established  from  widely  observed  proper¬ 
ties  of  turbulent  flows.  Their  values  are: 

3  =  3/20  B*  =  9/100 

a  =  1/2  a*  =  1/2 

Y#  =  {  1-(1-X2 )exp  (-ReT)f  (  (11) 

YY*  =  10/9  |  1-(1-X2)exp  (-ReT/2)f 
X  =  1/11 
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and  the  quantity  ReT  is  the  Reynolds  number  of  the  turbulence  defined 
in  terms  of  the  dependent  variables  by: 

ReT  =  pe/(ojy)  (12) 

Before  proceeding  to  a  discussion  of  boundary  conditions  in  the  next 
subsection,  a  key  point  about  the  equations  of  motion  is  worthy  of 
note.  In  the  mean-energy  equati->n  (4),  note  that  the  last  term  on 
the  right-hand  side  generally  is  omitted  by  turbulence  modelers.  The 
origin  of  this  term  is  obvious  when  the  time-averaged  equation  for 
total  energy,  E,  as  defined  in  Equation  (5)  is  derived.  If  this  term 
is  omitted,  there  is  no  guarentee  that  total  energy,  E,  will  be  con¬ 
served  in  the  flowfield,  even  when  e  <<  CvT.  This  is  true  because 
its  omission  can  give  rise  to  nonphysical  numerical  oscillations  whose 
energy  source  is  unrelated  to  the  physical  energy  of  the  flowfield. 

2.2  BOUNDARY  CONDITIONS 

As  noted  in  the  Introduction,  in  our  computation  of  the  flowfields 
of  interest,  we  Integrate  the  equations  of  motion  through  the  viscous 
sublayer.  Thus,  we  must  specify  boundary  conditions  for  the  various 
dependent  variables  appropriate  to  a  solid  boundary.  Denoting  distance 
normal  to  the  surface  by  y,  we  impose  the  no-slip  velocity  boundary 
condition  at  y=0.  Additionally,  for  all  compressible  applications 
considered  in  this  study,  the  surface  is  adiabatic.  Thus,  for  the 
mean-flow  quantities  we  have: 


u^  =  0 


8T/9y  =  0 


y  =  0 


Turning  now  to  the  turbulence  propoerties,  an  additional  consequence 

of  the  no-slip  boundary  condition  is  that  the  turbulent  mixing  energy 

(which  is  proportional  to  the  square  of  fluctuating  velocity  components) 

2 

vanishes  at  a  solid  boundary.  Also,  as  noted  by  Wilcox  and  Rubesin  , 
the  turbulent  dissipation  rate  approaches  a  known  analytical  form  in 
the  sublayer.  For  a  perfectly  smooth  surface  we  have: 
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where  subscript  w  denotes  surface  value. 


Equation  (15)  illustrates  a  key  point  about  the  turbulence-model 
equations  which  is  worthy  of  note.  Remembering  that  the  characteristic 
length  scale  of  the  turbulence,  i,  is  inversely  proportional  to  w 
approaching  a  solid  boundary  implies  that  l  becomes  very  small.  This 
is  consistent  with  the  physics  of  turbulence  which  finds  the  constraint 
of  a  solid  boundary  admitting  only  the  smallest  eddies.  Thus,  even 
with  the  great  advantages  attending  the  use  of  long-time  averaging 
(compared  to,  say,  a  large-eddy  simulation),  we  still  must  deal  with 
the  problem  of  resolving  very  small  eddies  in  a  finite-difference  com¬ 
putation.  This  "conservation  of  difficulty"  stands  at  the  center  of 
the  numerical  woes  attending  solution  of  the  turbulence-model  equations. 


Finally,  as  will  be  illustrated  in  the  next  subsection,  the  dissipation 

_p 

rate  is  singular  (although  less  strongly  than  y  )  far  above  the  vis¬ 
cous  sublayer.  Thus,  the  numerical  difficulties  we  have  encountered 
are  not  caused  entirely  by  our  treatment  of  the  sublayer.  Rather,  all 
advanced  turbulence  models  share  the  same  unpleasant  (from  a  numerical 
standpoint)  behavior. 


2.3  TRUNCATION  ERRORS 

Virtually  all  finite-difference  numerical  schemes  used  for  solving 

the  time-averaged  Navier  Stokes  equations  employ  central  differences. 

Usually,  such  schemes  are  of  second-order  accuracy  in  both  time  and 

space.  In  our  work,  for  example,  we  have  used  variants  of  the  MacCormack 

4-6 

time-splitting  method  which  uses  second-order  accurate,  central 
differences.  Even  for  relatively  crude  meshes,  such  a  scheme  is  sat¬ 
isfactory  when  the  various  dependent  variables  are  analytic  through¬ 
out  the  flowfield.  However,  when  one  of  the  dependent  variables  is 
nonanalytic,  such  as  oj,  truncation  errors  can  become  very  large. 

For  example,  Figure  1  compares  the  computed  near-surface  behavior  of 
w  with  the  theoretical  limiting  form  given  in  Equation  (15) •  The  profile 


« 


6 


Figure  1.  Comparison  of  computed  near— wall  dissipation— rate  profile  with  the 
asymptotic  limiting  behavior j  Mach  3  flat -plate  boundary  layer. 
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shown  is  for  a  Mach  3,  adiabatic -wall  boundary  layer.  Long  experience 
with  the  model  equations  has  shown  that,  for  a  wide  range  of  flows 
both  incompressible  and  compressible,  very  close  to  the  surface 
(y+  =  u  y/v  <  4)  Equation  (15)  is  valid.  As  shown,  the  computed 
values  of  w  are  20%  -  30%  higher  than  the  theoretical  values.  Our 
experience  with  these  model  equations  has  shown  that  this  large  an 
error  in  u>  has  a  large  effect  on  predicted  flow  properties  throughout 
the  boundary  layer. 


The  reason  such  truncation  errors  persist  far  above  the  sublayer  be¬ 
comes  obvious  upon  inspection  of  the  limiting  form  of  w  in  the  loga- 

*7 

rithmic  (law-of-the-wall )  region.  As  shown  by  Saffman  and  Wilcox' , 
approaching  the  sublayer  from  above  we  have 


ut  *yp 

/$*  Ky 


as 


(16) 


where  k  =  0.4l  is  Karman's  constant  and  u^  is  friction  velocity.  Thus, 
even  far  above  the  viscous  sublayer  w  is  nonanalytic  and  the  use  of 
central  differences  is  attended  by  nontrivial  truncation  errors. 


The  source  of  the  numerical  inaccuracy  becomes  obvious  if  we  consider 
the  magnitude  of  the  errors  introduced  when  using  central  differences 

p 

to  compute  derivatives  of  u)  close  to  the  surface.  Because  we  solve 
the  equations  of  motion  in  conservation  form,  the  finite-difference 
formulation  actually  computes  ptoz  which  (ignoring  density  for  simplicity) 
exhibits  strongly  singular  behavior,  viz, 
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It  is  easy  to  show  if  we  use  the  exact  values  of  u>2  at  two  adjacent 
mesh  points  y^  and  yj+1>  central  differences  yield: 
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Because  of  the  fine  resolution  close  to  the  surface  (for  the  point 
nearest  the  surface  the  computation  of  Figure  1  has  y+  =  H)  the 
smallest  value  of  e  for  the  mesh  used  is  1/5*  Using  this  value  for 
e  in  Equation  ( 1 8 )  yields 


(— ) 
v3y  ; 
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exact 
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Thus,  using  conventional  central  differences  yields  more  than  a  2555 
error  when  computing  derivatives. 

In  order  to  illustrate  the  point  that  the  truncation  errors  persist 
far  beyond  the  sublayer.  Figure  2  shows  results  of  a  computed  w 
profile  obtained  from  solution  of  a  simplified  version  of  Equation 
(7).  Specifically,  we  have  dropped  all  but  the  dissipation  and 

molecular  diffusion  terms  so  that,  for  incompressible  flow,  the 

equation  simplifies  to 

v  =  gw3  (21) 

dy2 


the  exact  solution  of  which  is 


0)  = 


20v 

By*" 


(22) 


The  figure  shows  results  of  two  computations  both  of  which  use  a  mesh 
similar  to  that  of  the  Navier-Stokes  computation  discussed  above. 

That  is,  mesh  points  normal  to  the  surface  are  spaced  in  a  geometric 
progression  with  progression  ratio  k  =  1.25.  In  the  first  computation 
the  value  of  w  at  the  mesh  point  nearest  the  surface  is  obtained  from 
Equation  (22).  As  shown,  the  local  error  close  to  the  surface  is  in 
excess  of  50$.  In  the  second  computation.  Equation  (22)  is  used  to 
prescribe  the  value  of  ui  for  the  10  mesh  points  nearest  the  surface 
(a  procedure  used  in  most  of  our  past  computations  provided  mesh  point 
number  10  lies  below  y+  =  4).  While  the  local  error  is  reduced  signif 
icantly  near  the  surface,  there  is  little  improvement  in  accuracy 
farther  from  y  =  0. 
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2.4  STABILITY  CONSIDERATIONS 


A  second,  less  obvious,  numerical  difficulty  becomes  evident  upon 
considering  the  stability  of  a  given  numerical  scheme  when  applied 
to  the  model  equations.  A  simplified  version  of  the  equations  for 
e  and  w2  which  illustrates  potential  problems  is  as  follows. 


3e 

at 


+  u 


3e 

3x 


=  Pe  + 


(23) 


In  Equation  (23),  ?  denotes  net  production,  viz,  production  minus 
dissipation.  In  the  case  where  P  is  constant,  it  is  easy  to  show 
that  even  for  a  typical  Implicit  second-order  accurate  time-marching 
scheme,  positive  P  leads  to  unconditional  instability'1'  in  the  classical 
sense.  This  is  unsurprising  however  as  solutions  to  Equation  (23), 
subject  to  certain  initial  and  boundary  conditions,  grow  exponentially 
in  time.  Thus,  even  though  local  errors  might  be  small  compared  to 
the  solution  itself,  they  will  also  grow  exponentially.  A  numerical 
analysis  of  Equation  (23)  performed  in  this  study  verified  this  in¬ 
stability  when  P  is  positive. 

A  key  shortcoming  of  the  model  equation  (23)  is  the  independence  of  e 
from  the  mean  kinetic  energy.  In  any  physical  flow,  there  will  not 
be  an  infinite  source  of  energy  as  implied  in  Equation  (23)-  Rather, 
the  coefficient  P  will  increase  at  the  expense  of  the  mean  kinetic 
energy  and  u  will  decrease.  To  simulate  this  situation  we  studied 
a  model  in  which  we  have 

P  =  pu2  - 

hw2  +  e  =  constant 

where  p  is  a  constant  and  an  equation  similar  to  Equation  (23)  is 

,  2j 

written  for  w  .  Using  MacCormack’s  explicit  scheme  ,  this  system 
displays  no  numerical  instability  for  Courant  numbers  less  than  one. 
Thus,  as  noted  in  Subsection  2.1,  if  we  are  not  careful  enough  in 
writing  our  turbulence  model  equations  to  insure  conservation  of 
total  energy,  unusual  numerical  instability  may  arise. 


(24) 

(25) 


However,  it  is  still  possible  for  the  production  term  to  cause  large 
numerical  excursions  from  the  steady  state,  especially  if  the  initial 
flowfield  is  quite  different  from  the  long-time  solution.  In  the  w 
equation,  for  example,  very  close  to  the  surface  we  should  find  a 
balance  between  dissipation  and  molecular  diffusion.  By  contrast  we 
have  often  observed  a  strong  (nonphysical)  coupling  between  the  un¬ 
steady  and  the  dissipation  terms  which  causes  u>2  to  become  negative. 
Hence,  the  presence  of  the  source  terms  can  present  numerical  problems 
quite  unlike  those  encountered  with  conventional  conservation 
equations. 


3-  A  NEW  ALGORITHM 
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In  this  section  we  present  a  new  algorithm  which  simultaneously 
eliminates  truncation  errors  and  effectively  removes  nonphysical 
coupling  between  source  and  unsteady  terms.  First  we  motivate  and 
present  the  algorithm;  then  we  apply  it  to  Equation  (21). 

3.1  MOTIVATION 

In  a  previous  study,  Wilcox1  tried  removing  the  near-wall  singularity 
analytically  by  making  the  straightforward  change  of  variables: 

w2  =  (Jj/y2*  (26) 

The  function  ^  is  analytic  as  y  -*■  0  and  presumed  a  more  appropriate 
quantity  to  compute.  However,  using  the  MacCormack  hybrid  explicit- 
implicit  method  Wilcox  simultaneously  (a)  improved  near-wall 
accuracy  and  (b)  introduced  large  (bounded)  numerical  oscillations. 
While  this  change  of  variables  eliminates  truncation  errors,  it 
apparently  aggrevates  numerical  stability  via  attending  alteration 
of  the  source  terms. 


To  provide  motivation  for  the  new  algorithm  which  will  be  presented 
in  the  following  Subsection,  it  is  instructive  to  consider  the 
limiting  form  of  the  w  equation  upon  approaching  y  =  0.  It  is  easy 
to  show  that  the  equation  has  the  following  form  for  a  constsnt- 
pressure  boundary  layer: 


3V 

3y2 


+  a  IhlI 
y  3y 


-  fi-  u>2  * 

y2 


(27) 


where  precise  details  of  the  coefficients  F,  p  and  q  are  unimportant 
(see  Equations  3^-35  below).  Of  greatest  importance  is  their  behavior 
close  to  the  surface  where  we  find: 


Viscous  Sublayer  (y+  <  10) 


F  -*■  0 
q  0 

p  +20 

y"w2->-constant 


as  y  +  0, 


y+  <  10 


(28) 
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Wall  Layer  (10  <  y+  <  100) 

F  0 

q  -*■  1 

p  -*•  4 

y2w2->-constant 


as  y  -*•  0,  10  <  y+  <  100 


(29) 


Inspection  of  Equations  (28)  and  (29)  shows  that  the  singularity 
changes  its  strength  across  the  sublayer/wall-layer  interface.  It 
futher  suggests  that  a  more-appropriate  transformation  might  reflect 
this  fact. 


3 . 2  THE  ALGORITHM 

The  observations  in  the  preceding  Subsection  suggest  the  following 
change  of  variables: 


01 


2 


yx 


(30) 


where  x  is  assumed  to  be  a  function  of  p  and  q  and  also  of  time,  t. 
Substituting  Equation(30)  into  Equation  (7)  yields 


If  ♦  } 

fPyHogy)^.  ,  +/likiLl<’p£>  -  PV* 

'■p  +  opeJ  3t  x  (  p  + 


ope 


-  l}x  + 


y2s 

p  +  ope 


(31) 

(32) 


where 

3  -  -  {«  +  2°(%)!}p“  <«> 
There  are  three  key  points  about  the  transformed  equations  worthy  of 
note.  First,  while  Equation  (31)  contains  a  nonconservation-like 
term,  viz,  the  term  proportional  to  3\p/d y,  the  term  shouldn't  destroy 
numerical  stability.  The  new  term  will  behave  more  like  an  extra 
convective  flux  and  should  be  far  less  troublesome  provided  it  is 
differenced  in  a  manner  consistent  with  that  used  for  the  standard 
convective  terms. 


Second,  the  coefficients  in  the  equation  for  x  are  directly  related 


» 


14 


r  v.  ... ;  t 


•  • 


I 


I 


» 


» 


t 


t 


t 


> 


► 


to  the  p  and  q  appearing  in  Equation  (27).  In  fact,  we  have 

y{|^(y  +  ape)  -  pv| 

^  U  +  ape 

P  =  -  ^ - 

H  U  +  ape 

Hence,  in  the  steady  state  (3x/3t  =  0)  we  have 


X2 3  +  (q  -  l)x  -  P  =  0 


so  that 


X  =  *s(l  -  q)  -  Vtl  -  q) 1  +  "4p 


(34) 

(35) 

(36) 

(37) 


where  we  have  selected  the  root  of  Equation  (36)  which  yields  negative 
X*  Substituting  the  values  of  p  and  q  from  Equations  (28-29)  into 
Equation  (37)  shows  ? c-v. 


A 


in  the  sublayer 

in  the  wall  layer 


(38) 


Thus,  the  transformation  captures  the  proper  nature  of  the  singularity 
in  both  near-wall  regions. 


The  final  key  point  is  that  we  have  a  method  which  requires  no  specific 
alteration  of  the  particular  numerical  method  used.  Central  differ¬ 
ences  should  be  satisfactory  for  computing  as  it  is  an  analytic 
function  as  y  -*■  0. 

In  summary,  the  change  of  dependent  variables  proposed  in  Equation  (30) 
has  the  following  three  desirable  features: 

1.  It  replaces  the  Pu2  type  source  term  with  a 
more  familiar  convective-type  term; 

2.  It  reflects  the  varying  strength  of  the  w2 
singularity  In  the  near-wall  region  of  a 
boundary  layer; 

3.  It  places  no  unusual  constaint  on  the  type  of 
numerical  algorithm  used. 
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In  principal,  this  transformation  has  very  desirable  properties. 

To  provide  a  definitive  numerical  accuracy  test,  the  next  Subsection 
shows  results  of  its  application  to  Equation  (21)  for  which  the  exact 
solution  is  known,  i.e..  Equation  (22). 

3-3  APPLICATION  WITH  A  SIMPLIFIED  MODEL 

t  As  our  first  application  of  the  new  algorithm,  we  consider  the  limit¬ 

ing  near-wall  form  of  the  equation  for  to,  viz.  Equation  (21).  This 
equation  provides  a  definitive  accuracy  check  as  the  exact  solution 
is  known  and  is  given  in  Equation  (22).  The  new  MacCormack^  implicit 
I  method  is  used  so  that  we  introduce  an  unsteady  term  to  Equation  (21) 

and  seek  the  long-time  solution.  Thus,  the  actual  equation  being 
solved  is 


t 


« 


» 


» 


> 


3w2  _  32oj2  r  3 

It  "  V  JyT  ~ 

(39) 

Introducing  the  transformation  defined  in  Equation  (30) 

following  equations  for  ip  and  x: 

leads  to  the 

ii  _  „  lii  .  ixz  21 
at  "  v  3y*  y  ay 

(40) 

(S!l2£l,|X  .  X2  .  x  .  eyWv 

(41) 

where  w  is  understood  to  be  computed  from  \p  and  x  according  to  Equation 
(30).  In  the  steady  state,  the  exact  solution  to  Equations  (4o)  and 
(41)  is 


<Ji  =  2Ov/0  (42) 

X  =  -4  (43) 

These  equations  have  been  solved  using  a  grid  with  equally  spaced 
mesh  points.  Such  a  grid  which  would  yield  even  larger  truncation 
errors  than  shown  in  Figure  2  if  we  solved  Equation  (39)  without  the 
transformation  to  \|>  and  x*  The  grid  was  also  constructed  to  insure 
that  the  nondimensional  y  coordinate,  viz, 

y  =  Uy/v  (44) 
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where  U  is  a  characteristic  velocity,  will  always  be  greater  than 
unity.  Early  numerical  experimentation  showed  that  Equation  (41) 
is  more  easily  solved  if  the  logy  term  on  the  left-hand  side  re¬ 
tains  the  same  sign  throughout  the  mesh. 

The  initial  profiles  for  and  x  were  varied  about  the  exact  steady- 
state  values  given  in  Equations  (42)  and  (43). 

In  all  computations,  no  more  than  50  timesteps  were  required  to  reach 
steady-state.  Regardless  of  the  initial  conditions,  the  solution 
converged  to  within  six  significant  figures  of  the  exact  solution. 

Thus,  at  least  for  this  simplified  test  case,  the  transformation 
defined  in  Equation  (30)  yields  greatly  improved  numerical  accuracy. 

Because  the  exact  solution  to  Equations  (40-41)  is  *  constant  and 
X  =  constant,  this  simplified  example  is  too  simple  to  warrent  pro¬ 
claiming  the  new  method  as  proven.  Nevertheler s ,  this  test  is  well 
defined  and,  if  the  method  failed,  there  would  be  little  point  in 
pursuing  it  any  farther.  In  the  next  Section,  we  put  the  new  algorithm 
to  two  rigorous  tests. 


4.  TEST  CASES 


In  this  section,  we  present  results  of  two  computations  using  the 
MacCormack  implicit  time-splitting  scheme  in  conjunction  with  the 
new  algorithm.  The  full  equations  of  motion  are  solved  for  two 
compressible  boundary  layers  with  a  full  time-averaged  Navier-Stokes 
computer  program.  The  first  application  is  to  a  Mach  2  laminar 
boundary  layer  and  the  second  to  a  Mach  3  turbulent  boundary  layer. 

4.1  MACH  2  LAMINAR  BOUNDARY  LAYER 

Our  first  application  addresses  a  Mach  2  flat-plate  laminar  boundary 
layer.  This  case  provides  a  significant  test  of  the  new  algorithm 
for  two  reasons.  First,  the  near-wall  w  profile  defined  in  Equation 
(15)  is  known  to  hold  through  a  significant  portion  of  the  boundary 
layer  so  we  have  a  rigorous  check  on  accuracy.  Second,  previous 
computations  have  often  had  a  nonphysical  transition  to  turbulence 
even  at  very  low  Reynolds  numbers. 

The  finite  difference  mesh  consists  of  rectangular  cells  with  12 
mesh  points  in  the  streamwise  direction  and  16  points  normal  to  the 
surface.  Mesh  points  are  spaced  in  a  geometric  progression  with  a 
progression  ratio  k  =  1.20.  Computation  is  initiated  from  uniform 
flow  with  the  upstream  boundary  just  ahead  of  the  plate  leading  edge. 
The  downstream  boundary  lies  at  a  plate-length  Reynolds  number  of 
4-loV  The  computation  was  run  for  50  timesteps.  Steady  flow  con¬ 
ditions  appear  to  have  been  obtained  after  about  35  timesteps  with 
the  quantity  x  taking  longest  of  all  dependent  variables  to  approach 
steady  state.  Computing  time  for  this  case  is  six  hours  on  a  TRS-80 
microcomputer. 

Figure  3  compares  a  computed  w  profile  with  the  theoretical  near-wall 
solution.  As  indicated,  the  two  solutions  are  in  very  close  agree¬ 
ment  up  to  about  y+  *  uTy/v  =  8.  Above  this  point  the  near-wall 
solution  is  not  expected  to  be  valid  as  convective  terms  become  im¬ 
portant.  Thus,  we  can  only  check  our  solution  in  the  asymptotic 


LAMINAR  BOUNDARY  LAYER 


Figure  3-  Comparison  of  computed  near-wall  dissipation-rate  profile  with 
the  asymptotic  limiting  behavior;  Mach  2  Laminar  boundary  layer 


limit  y  -*•  0.  However,  this  is  precisely  the  limit  of  greatest 
concern  and  the  overall  agreement  is  excellent. 

At  no  point  in  the  computed  flowfield  is  there  any  evidence  of 
turbulent  energy  amplification  which  might  lead  to  nonphysical 
transition.  The  new  algorithm  thus  appears  to  eliminate  this 
undesired  computational  result. 

Because  the  boundary  layer  remains  laminar  in  this  computation, 
the  two  turbulence  model  equations  play  a  passive  role.  To  provide 
an  even  more  rigorous  test  we  must  now  turn  to  a  case  in  which  the 
model  equations  play  an  active  role.  In  our  next  application  we 
address  such  a  case. 

4.2  MACH  3  TURBULENT  BOUNDARY  LAYER 

In  this  application  we  consider  a  Mach  3  turbulent  boundary  layer. 
Again,  the  mesh  consists  of  12  points  in  the  streamwise  direction 
and  16  points  normal  to  the  surface.  The  mesh  point  nearest  the 
surface  lies  at  y+  =  2  and  the  grid  points  are  placed  in  a  geometric 
progression  with  k  =  1.25.  Initial  profiles  are  obtained  from  a 
boundary-layer  program. 

The  length  of  the  mesh  is  22  boundary-layer  thicknesses  and  Reynolds 

6  6 

number  based  on  plate  length  extends  from  1.8-10  to  2.5-10  .  The 
momentum-thickness  Reynolds  number  is  2000  at  the  upstream  boundary. 

The  main  shortcoming  of  previous  Navier-Stokes  solutions  of  the 
model  equations  is  that,  in  addition  to  w  being  20$  -  30$  too  large 
in  the  sublayer,  skin  friction  generally  Is  15$  low.  Hence,  in 
addition  to  checking  the  accuracy  of  the  w  profile,  we  must  also 
check  the  skin  friction. 

Figure  4  compares  computed  and  theoretical  a  profiles  for  the  near¬ 
wall  solution  up  to  y+  =  100.  As  shown,  the  numerical  solution  is 
within  8$  of  the  theoretical  solution  below  y+  =  10.  Between  y+  =  10 
and  100,  the  numerical  solution  asymptotically  approaches  the  wall- 
layer  profile  (Equation  16).  Thus,  u>  has  been  computed  far  more 
accurately  than  in  previous  Navier-Stokes  calculations  (cf  Figure  1). 
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Figure  4.  Comparison  of  com*  Jted  near-wall  dissipation-rate  profile  with  the 
asymptotic  limiting  behavior;  Mach  3  turbulent  boundary  layer. 


The  skin  friction  coefficient  is  about  5%  lower  than  in  the  correspond¬ 
ing  boundary-layer  computation.  This  is  somewhat  surprising  as  the 
w  profile  is  actually  underpredicted  near  the  surface.  Our  earlier 
arguments  would  predict  a  skin  friction  coefficient  which  is  some¬ 
what  overpredicted.  The  coarseness  of  the  mesh  may  be  responsible. 
Nevertheless,  this  application  indicates  a  great  improvement  in 
numerical  accuracy. 

This  computation  required  125  timesteps  to  reach  steady  flow  conditions, 
partly  because  the  timestep  was  smaller  than  in  the  laminar  case  and 
partly  because  x  converged  to  its  steady-state  value  very  slowly. 

The  overall  computing  time  was  15  hours  on  a  TRS-80  microcomputer. 
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5.  RESULTS  AND  CONCLUSIONS 

I  The  new  algorithm  appears  to  greatly  improve  numerical  accuracy  of 

a  time-marching  method  when  used  to  compute  a  quantity  such  as  uj 
which  behaves  nonalytically  approaching  a  solid  boundary.  This 
is  particularly  true  when  a  very  crude  finite-difference  mesh  is 
►  employed.  Although  the  algorithm  has  been  tested  for  relatively 

simple  flows,  we  have  no  reason  to  suspect  it  will  fail  in  more 
complex  applications. 

^  An  especially  desirable  feature  of  the  new  algorithm  is  its  indepen¬ 

dence  of  the  actual  numerical  method  being  used.  In  essence,  we 
have  found  an  analytical  method  for  removing  a  singularity  of  vary- 

_ii  -2 

ing  strength  (y  very  close  to  the  surface,  y  farther  from  the 
t  surface).  In  so  doing,  the  actual  equations  which  are  solved 

numerically  can  be  accurately  solved  by  any  conventional  method 
using  central  differences. 

^  Results  presented  in  Sections  3  and  4  warrent  further  testing  and 

development  of  the  technique.  The  most  immediate  tests  should 
be  for  more  complex  flows  such  as  shock  induced  turbulent  boundary- 
layer  separation.  More  development  effort  is  needed  to  generalize 
^  the  transformation  defined  in  Equation  (30)  for  arbitrary  geometry. 

After  this  testing  and  development  has  been  done  we  should  be  ready 
to  proceed  toward  our  original  goal,  viz,  that  of  generating  accur¬ 
ate  numerical  solutions  for  general  flow  applications  such  as  those 
i  occurring  in  axial  turbomachines. 
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CONTRACT  OVERVIEW 


The  overall  object  of  this  project  has  been  to  develop  a  method  for 

computing  flow  through  the  rotor  passages  of  highly-loaded  axial 

turbomachines.  Initially  the  technical  approach  was  to  use  the  MacCormack 
5 

hybrid  method  to  compute  the  inviscid  flow  and  to  use  a  three-dimen¬ 
sional  boundary-layer  program  to  resolve  the  boundary  layers  up  to 
separation.  As  the  project  progressed,  it  became  evident  that  coupling 
the  inviscid  flow  program  with  the  boundary-layer  program  was  proving 
far  more  difficult  than  anticipated  because  of  numerical  problems 
associated  with  (a)  obtaining  smooth  pressure  gradient  input  for  the 
boundary-layer  program  and  (b)  numerically  accurate  boundary-layer 
profiles  without  using  inordinate  numbers  of  mesh  points  normal  to 
blade  surfaces.  Noting  the  great  promise  offered  by  the  new  MacCormack 
implicit  method^,  we  decided  to  eliminate  the  coupled  inviscid/ 
boundary-layer  computational  procedure  in  favor  of  a  full  time-averaged 
Navier-Stokes  computation. 

A  sophisticated  mesh  generation  procedure  was  developed  for  rotor 
passages  and  all  program  modifications  were  made  to  incorporate  viscous 
stresses  in  the  program. 

At  this  point  we  found  that  serious  truncation  error  and  numerical 
stability  problems  were  present  in  the  computer  program  which  pre¬ 
cluded  definitive  testing  of  the  method.  These  problems  are  delineated 
in  Section  2  of  this  report.  In  order  to  analyze  and  eliminate  these 
problems,  the  study  reported  herein  was  conducted.  Although  the 
algorithm  devised  in  this  study  has  not  been  incorporated  in  our  turbo¬ 
machine  program,  we  have  little  doubt  that  it  will  resolve  the  dif¬ 
ficulties  encountered. 

An  additional  task  was  added  to  the  contract  to  perform  computations 
for  the  1981  Stanford  Conference  on  Complex  Turbulent  Flows.  Results 
of  those  computations  appear  in  Reference  4  of  the  Bibliography. 
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